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Abstract 

Since 1996, a hybrid experiment consisting of the emulsion chamber and 

burst detector array and the Tibet-II air-shower array has been operated at 

Yangbajing (4300 m above sea level, 606 g/cm^) in Tibet. This experiment 

can detect air-shower cores, called as burst events, accompanied by air showers 

in excess of about 100 TeV. We observed about 4300 burst events accompanied 

by air showers during 690 days of operation and selected 820 proton-induced 

events with its primary energy above 200 TeV using a neural network method. 

Using this data set, we obtained the energy spectrum of primary protons in 

the energy range from 200 to 1000 TeV. The differential energy spectrum 

obtained in this energy region can be fitted by a power law with the index 

of -2.97 lb 0.06, which is steeper than that obtained by direct measurements 

at lower energies. We also obtained the energy spectrum of helium nuclei at 

particle energies around 1000 TeV. 

PACS numbers : 98.70Sa, 95.85Ry, 96.40De, 96.40Pq 
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I. INTRODUCTION 



Shock acceleration at supernova blast waves gives a good explanation of the origin of 
the bulk of cosmic rays. It may be well accepted that cosmic rays below about 10 TeV are 
predominantly due to the explosion of stars (supernova explosion) into the normal interstellar 
medium, while particle acceleration at supernova remnants (SNR's) has an upper limit of 
about 100 TeV 0,0] • Also, there is an argument that the cosmic rays from near 10 TeV 
to several times 1000 TeV very likely originate in the explosion of massive stars into their 
former stellar wind ||^. These processes have been examined to be able to explain the 
cosmic ray spectra fairly well up to the highest energy where abundances are known . For 
energies beyond about 1000 TeV, however, there is no consensus. On the other hand, ground- 
based air-shower experiments observe cosmic rays with energies up to ~100 EeV (10^° eV). 
Measurements so far reported suggest that the slope of the all-particle spectrum in the 
energy range of about 100 - 1000 TeV is somewhat flatter than that observed at lower 
energies, while at higher energies over several times 1000 TeV the energy spectrum becomes 
steeper with the slope of about -3.0. The break in the overall spectrum at around 1000 TeV 
is often referred to as the "knee" in the spectrum. 

Clearly, the knee of the primary cosmic ray spectrum has its origin in the accelera- 
tion and propagation of high-energy cosmic rays in our Galaxy. The acceleration model 
by supernova blast waves leads to the formation of a power-law spectrum of particle en- 
ergies with the index of about -2 at sources [Q, and plausible propagation models of their 
confinement by galactic magnetic fields and of their eventual escape from our Galaxy can 
explain well a steeper power-law spectrum than that at the source region [^], suggesting a 
rigidity-dependent bending for different cosmic ray composition. Within the framework of 
this picture the average mass of primary cosmic rays before the knee should increase with 
increasing primary energy. In other words, the knee composition becomes heavy dominant 
as the proton spectrum may first bend at an energy of about 100 TeV, corresponding to a 
maximum energy gained by shock acceleration at SNR's |]1|. 
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While the origin of cosmic rays with energies beyond the knee is still in controversy, 
observations of cosmic rays in such a high energy region may naturally stand in need of 
other acceleration mechanisms [^|^ or new cosmic ray sources 0,^. Among those, one of 
the most promising models may be that the cosmic rays come from extra-galactic sources 
such as active galactic nuclei [Q, though the evidence is far from convincing. However, such 
an extra-galactic source model should predict proton-enriched primary composition around 
and beyond the knee. 

Thus, measurements of the primary cosmic rays around the knee are very important and 
its composition is fundamentally input for understanding the particle acceleration mech- 
anism that pushes cosmic rays to very high energies. Among various primary particles, 
protons hold the key to the situation and its spectrum provides major constraints on the 
model parameters of the origin of high-energy cosmic rays. Because of extremely low and 
steeply decreasing flux at high energies, however, direct measurements of primary proton 
spectrum on board balloons are still limited in the energy region lower than a few hun- 
dred TeV. In a recent report by the JACEE group |]Tn| it was concluded that the proton 



spectrum as well as the helium spectrum are consistent with power laws with no spectral 
breaks, meaning that there is no bending up to the highest energy they measured (about 
800 TeV). However, this is a surmise based on statistically sparse data, so more studies 
are required. On the other hand, most studies on the cosmic ray composition around the 
knee have been carried out with ground-based instruments that can observe the various air- 
shower parameters. Recently, for example, measurements of muon content in each air-shower 
TT| or muons in the deep underground |[T^JT3[] , measurements of the lateral distribution of 



air shower Cherenkov lights or the maximum depth of shower development using air 



Cherenkov telescopes and multiparameter measurements of air showers [|T^ have been 
carried out. However, the results obtained by these methods have been derived by indirect 
ways that may strongly rely on how the observed quantities depend on the composition, on 
the precision of the measurements, and on the air-shower and detector simulations as well. 
Therefore, the conclusions sometimes differ with experiments considerably. 



Within the ground-based experiments those which set up at higher altitudes are prefer- 
able. The reasons include, first, the observation level is close to the maximum of the shower 
developments induced by cosmic rays with energies around the knee, so that the energy 



determination is more precise and less dependent upon the unknown composition [|T^] ; sec- 
ond, the higher energy flux in the core region of air showers can be observed with emulsion 
chambers or burst detectors ||l8l . High-energy air-shower cores are sensitive to the intensity 
of protons in the primary cosmic rays and also to the composition around the knee. 

A hybrid-experiment of emulsion chamber and air-shower array at high altitude has a 
great advantage for studying the composition of primary particles at the knee energy region 
T9| , |20[| . In a previous paper we have developed a method to study the primary cosmic 



ray composition with a hybrid detector of the emulsion chamber and air-shower array based 
on a Monte Carlo simulation. It is shown there that an artificial neural network (ANN) 
can be used as a classifier for the species of primary particles since high-energy air-shower 
cores accompanying air showers are characterized by several parameters and that such a 
hybrid experiment is powerful enough to select the events induced by protons in the knee 
energy region. We have applied a three layered feed forward neural network with a back- 
propagation learning algorithm to the data obtained with the Tibet burst detector and the 



air-shower array [P2 



Here, we report our study on the primary proton spectrum using the data obtained with 
the Tibet burst detector and air-shower array. The experiment, including the apparatus, its 
performance and data selection, is described in Sec. II. Air-shower simulations to compare 
with the experimental data are described in Sec. HI. The ANN used is briefly introduced in 
Sec. IV. Section V is devoted to the results and discussions and a brief summary is given in 
Sec. VI. 



II. EXPERIMENT 
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A. Apparatus 



We started a hybrid experiment of the emulsion chamber, the burst detector and the 
air-shower array (Tibet-II) at Yangbajing (4300 m above sea level, 606 g/cm^), Tibet in 
1996 IQ. The Tibet-II array consists of 221 scintillation counters of 0.5 m^ each of which 
are placed on a 15 m square grid, and which has been operated since 1995. Any fourfold 
coincidence in the detectors is used as the trigger condition for air-shower events. Under 
this condition the trigger rate is about 200 Hz with a dead time of about 12% for data 
taking. The energy threshold is estimated to be about 7 TeV for proton-induced showers. 
The precision of the shower direction determination is about 1°, which has been confirmed 



by observing the Moon's shadow . The main aim of Tibet-II is to search for gamma ray 
point sources at energies around 10 TeV. But it can also be used for the measurement of 
the all-particle spectrum of cosmic rays 0, and for the study of topics in the knee region by 
providing information on the shower size, direction, core position, and arrival time of each 



air-shower event to the core detectors 22.24 



The emulsion chambers and the burst detectors are used to detect high-energy air-shower 
cores accompanied by air showers induced by primary cosmic rays with energies above 10^^ 
eV. They are separately set up in two rooms as shown in Fig. ^ and placed near the center of 
the Tibet-II array. A basic structure of each emulsion chamber used here is a multilayered 
sandwich of lead plates and photosensitive x-ray films [|l^. Photosensitive layers are set 
every 2 cascade units (c.u.) (here, 1 c.u. is taken to be 0.5 cm) of lead in the chamber as 
shown in Fig. 0. There are 400 units of emulsion chamber, each with an area of 40 cm x 50 
cm with the total thickness of 15 c.u., giving the total sensitive area of 80 m^, and 100 units 
of burst detectors each with an effective area of 160 cm x 50 cm. Four units of the emulsion 
chamber are set above one unit of the burst detector. A 1 cm iron plate is set between the 
emulsion chambers and burst detectors. 

Each burst detector consists of a plastic scintillator with the size of 160cm x 50cm 
and thickness of 2 cm, and four photodiodes (PD's) are attached at four corners of each 
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scintillator to read light signals generated by shower particles produced in the lead and iron 
absorber above the detector. Using the analog-to-digital converter (ADC) values from four 
PD's the total number (i.e., burst size, A*";,, for each burst detector) and the position of the 
number-weighted center of all shower particles that hit a burst detector can be estimated. 
The response of the burst detector is calibrated using electron beams from an accelerator 
and cosmic ray muons. The performance of the burst detector and the calibration using the 
electron beams are briefly summarized in Appendix A. It is confirmed that the burst size 
capable of measuring with each detector ranges from 10^ to 3 X 10^ roughly corresponding 
to showers with energies ranging from ~ 2 to ~ 300 TeV. 

A burst event is triggered when any twofold coincidence of signals from four PD's of a 
burst detector appears. Using the burst detector array shown in Fig. |l], the electromagnetic 
components in the air-shower cores can be measured in the area within a radius of several 
meters. The coincidence of a burst event and an air-shower event is made by their arrival 
times, and the coincidence of a burst event and a family event observed in the emulsion 
chamber is made by their positions and directions (A burst event and its accompanying 
air-shower have the same direction.). 

In the following analysis we use only the data obtained from all burst detectors and 
the Tibet-II array, while the emulsion chamber data will be reported elsewhere in the near 
future. 



B. Data analysis 

The data set of the burst events analyzed in this paper was obtained during the period 



from October 1996 through June 1999 p4|. First we scan the target maps of all events by 
the naked eye. Some events showing a systematic noise configuration were ruled out during 
the first scanning. An example of the burst detector event is shown in Fig. ^ where the 
size of the rhombus is logarithmically proportional to the burst size. A remarkable lateral 
distribution in the event pattern is seen. 
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Here, for convenience we introduce a "TOP detector" for each burst, which is defined 
as a detector recording the highest burst size among all fired burst detectors. Furthermore, 
since all the burst detectors are separately set up in two sections with a fairly large distance 
of 9 m as shown in Fig. ^ we call the section containing the TOP detector the "TOP 
section," and the other the "OTHER section" 

We first examined whether the burst detectors located far from the TOP detector still 
contain signals. For this, we divided the all burst events into five groups according to the 
case that the TOP detector in each event is in the first, second, third, fourth, and fifth 
column and then the size distribution observed with burst detectors in the OTHER section 
was obtained for each group. If almost all bursts observed in the OTHER section are signals, 
their size distribution must be different from event to event because they have different core 
distances. However all these five curves are of the same distribution, as seen in Fig. |^. 

This may strongly suggest that the bursts recorded in the detectors far away from the 
TOP detector, i.e., air-shower core, by more than 10 m are mostly formed by some noises, 
and its (equivalent) burst size (Nf,), which is estimated from the ADC value, is always smaller 
than 3 x 10^ under our experimental conditions as seen in Fig. ^. Here, the burst size of 
3 X 10^ corresponds to a few to 10 TeV for a single gamma ray or a single electron incident on 
the surface of the emulsion chamber. These noises may be mostly induced by an incomplete 
ground connection of the detectors to the earth. Hence, we subtracted the background in the 
TOP section assuming that the same background as in the OTHER section should appear 
in the TOP section and they randomly distribute in position. 

After the background subtraction, for a further analysis we made the data set by imposing 
the following conditions on the observed events : (1) Size of a TOP detector, nI°^ > 10^ ; 
(2) size of any non TOP detector, jY^"°'^~*°P > xo^ ; and (3) number of fired detectors with 
Nb > 10^ Nbd > 1- The total burst size for each burst event is defined as Y^N^,, where the 
summation is over all fired detectors with A^^ > 10^. 

5627 events are selected by these criteria, and among them 4274 events are accompanied 
by air showers with A'e(shower size) > 10^'^, which are recorded by the Tibet-II array. 
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The time intervals between two neighboring events are analyzed, and a good exponential 
distribution is seen, indicating a good randomness of this data sample. The effective running 
time of this experiment was estimated to be 689.5 days. Since the burst detector array was 
triggered separately with the Tibet-II array that has a 12% dead time, this value is taken 
into account when we calculate the intensity and the number of effective events. 



III. MONTE CARLO SIMULATION 

An extensive Monte Carlo simulation was carried out to simulate the cascade devel- 
opments (air showers) of incident cosmic rays in the atmosphere and the burst detector 
responses. To generate air-shower events in the atmosphere, we used two simulation codes, 
CORSIKA (+QGSJET interaction model) H and COSMOS pB], both of which are widely 



used in air-shower experiments. We also used an EPICS code ||2^ to simulate electromag- 
netic cascade showers in the detector. In this simulation, the detector performance, trigger 
efficiency of detectors, and effective area are adequately taken into account, based on the 
experimental data. 



A. Primary composition 

Primary particles we assumed were classified into seven species as proton (abbreviated to 
P and mass number=l), helium (He, 4), light nuclei (L, 8), medium nuclei (M or CNO, 14), 
heavy nuclei {H, 25), very heavy nuclei {VH, 35), and iron group (Fe, 56). The absolute flux 
of each composition was fitted to that obtained by direct measurements in the energy region 
around 1 ~ 10 TeV. The extrapolation to higher energies depends on the slopes of energy 
spectra and their bending points. As in our previous studies |^T],^, the heavy dominant 
(HD) and proton dominant (PD) models were examined. In HD (PD) the power indices were 
assumed to be 2.75 (2.65) for P, 2.65 (2.65) for He, 2.70 (2.70) for L, 2.52 (2.60) for M, 2.60 
(2.60) for H, VH, and 2.4 (2.60) for Fe, respectively. The bending energy was assumed to 
be proportional to the charge number and for protons to be 100 TeV in HD, while 2000 TeV 
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for all compositions in PD. The fractions of the proton component to the total at 100 and 
1000 TeV are 23 and 11 % in HD, and 40 and 39 % in PD, respectively. In both cases, the 
absolute intensity of all particle spectrum was normalized so as to be able to reproduce the 



Tibet and other experimental data well EM. The energy spectra of respective components 



assumed in the HD and PD models are summarized in Appendix B. 

B. Simulation procedure and simulation data 

Primary particles at the top of the atmosphere were sampled isotropically for the zenith 
angles within 45°. The minimum sampled energy of primary protons was set to 79 TeV and 
for other nuclei their minimum energies are determined so as to keep their contributions 
from lower energies to be less than 1%. All shower particles were followed till 5 GeV by a 



full Monte Carlo method and then till 1 GeV by the thinning method p5| , |28[| . The shower 
particles lower than 1 GeV were found to give minor contribution to the burst size since 
they are absorbed in the lead and iron. The air-shower size of each event was obtained using 
the data calculated by the thinning method. 

Each air-shower core which contains all shower particles with energies above 1 GeV was 
thrown on the burst detector array. Cascade developments of these shower particles in the 
burst detectors were calculated by use of the analytical formula which can well fit the full 
Monte Carlo simulation data obtained by EPICS ||2^. The selection of simulated burst 
events and their analysis were done under the same conditions as used for the experiment. 

The events were selected from the simulation data by imposing the same criteria as the 
experiment, and we obtained 4 xlO^ events (9200) for the CORSIKA+HD model (COS- 
MOS+HD). Among those selected events, 50% (48%) were induced by protons, 19% (17%) 
by helium, 17% (15 %) by L — CNO, and 14% (20 %) by other heavy nuclei, respectively, 
while for CORSIKA+PD, 2 xlO^ events were obtained and the primary ratios are 74%, 
16%, 7.5%, and 2.5%, respectively. The number of simulated events are 15 times as many 
as the experimental data. It may be worth noting here that the proton-induced events are 
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preferentially selected when air showers are tagged by high-energy cores. That is, even if the 
primary is heavy-enriched, almost half of the observed events selected by the above criteria 
are induced by protons. This is the reason why we can obtain the primary proton flux from 
this experiment successfully. 

Each event obtained can be characterized by the following three parameters : (1) Total 
burst size, J2 (2) total number of fired burst detectors, A^'^^i ; and (3) shower size, Ng. 
Among the three parameters, J2 and A^e are fairly sensitive to the primary composition, 



as discussed in the previous paper p4[. The scatter plots between J^Nb and for the 
CORSIKA+HD model are shown in Fig. ^. It is seen that the events with smaller Ne 
and larger J2 are mostly generated by protons. We use these simulation events in the 
following analysis. 

IV. ARTIFICIAL NEURAL NETWORK ANALYSIS 

As discussed above, the burst events accompanied by an air-shower are well characterized 
by the air-shower size, burst size, and the number of fired burst detectors. In this experiment, 
it is also noted that proton-induced events can be characterized by small air-shower size and 
large burst size, while those induced by heavy nuclei have the opposite character as their 
production height is relatively high in the atmosphere because of shorter mean free path 
than protons. Based on these facts, a simple multivariant analysis was introduced to select 



proton-induced events pO[. However, air-shower events are very complicated and it is not 
always obvious what data selection (or cuts) optimally enhance the signal (proton induced 
events) over the background. Neural networks may be an effective tool since they are ideal 
for separating patters into categories (e.g., signal and background). We can train a network 
to distinguish between signal and background using many parameters to describe each event. 
The network computes a single variable that ranges from zero to one and if the training is 
successful the network will output a number near zero for a signal event and near 1 for 
a background event. Hence, a single cut can be made on the network output which will 
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enhance the signal over the background. 

Usually, in a classification problem like the separation of proton-induced events and oth- 
ers, a set of p events with kmax observed variables each, described by the input vector {x^^^ 
= {xi,X2, ...,Xk^^^) has to be assigned to output categories yi using a set of classification 
functions t/i = Fi{{x}). For an example, a separation between signal and background events 
may be based on a one-dimensional output yi with the desired value for proton events and 
1 for other events. 

For a feed forward artificial neural network (ANN) with one layer of hidden units the 
following form of Fi is often chosen : 

Fiii^}) = aC^i^ijdC^'i^jkXk + Oj) + 9i), (4.1) 

j k 

which corresponds to the architecture of Fig. ^. Here, the weights Wij and wjk are the 
parameters to be fitted to the data distributions, and 6i and 6j are the thresholds which are 
generally omitted in the description as they can always be treated as weights 6i = Wio with 
Xo = 1. 

g{x) is the nonlinear neuron activation function, typically of the form (sigmoid function) 

^(x) = i[l + tanh(|;)], (4.2) 

where T is a parameter called temperature which is usually set to 1. 

The bottom layer (input) in Fig. |^ corresponds to sensor variables Xk and the top layer 
to the output features yi (the classification function Fi). The hidden layer enables nonlinear 
modeling of the sensor data. The great success of neural networks is mainly based on the 
derivation of an iterative learning algorithm based on gradient descent, the so-called back- 
propagation algorithm, and the weights Wij and wjk are determined by minimizing an error 
measure of fit, e.g., a mean-square error 

E = lT.iy^-t?r (4-3) 

p,i 

between y^ and the desired feature values tj with respect to the weights and (p) is an element 
of the training data sample. 
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Changing uoij by gradient descent corresponds to 

Aujij = -r]d,hj + aAujf/ (4.4) 
for the hidden to output layers, where 6i is given by 

5i = iVi - ti)g'{J2 ^ijhj)- (4-5) 
j 

Correspondingly, for the input to hidden layers one has 

Aujk = -vYl ^iAa'CZl ^ji^i)^k + OiAufjf. (4.6) 

i I 

In Eqs. (4) and (6) 1] is a. learning strength parameter which controls the speed of weight 
adjustment, and so-called momentum terms aAuj"^'^ and aAu^lf are included to damp out 
oscillation. A constant a determines the effect of the previous weight change. When no 
momentum terms are used, it takes a long time before the minimum has been reached with 
a low learning rate, whereas for high learning rates the minimum is never reached because of 
the oscillations. For a detailed description of the network technique, the back-propagation 
algorithm and modifications of the learning rule, see, e.g., |^ 



In this analysis, each data set is divided into two parts ; one that is used for training the 
network (training data set) and the other that is used for testing the ability of the network 
(test data set). Then, the whole training data sample is repeatedly presented to the network 
in a number of training cycles. After the network training an independent test data is used 
to check whether the network is able to generalize the classification to the data observed by 
our experiment. 

In this work we used a three-layered feed forward network as classifier of the species of 
primary particles. That is, this network contains three parameters as input neurons, ten 
hidden nodes, and one output unit and is abbreviated to a 3:10:1 network. Three parameters 
as input variables are ; (1) Air shower size Ng ; (2) the number of fired burst detectors Nbd, 
; and (3) sum of the size of fired burst detector J2 Nb. 

These are obtained for each event with the detector system consisting of the Tibet-II 
array and 100 burst detectors each with an effective area of 160 cm x 50 cm. The weights 
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in the network were initialized as uniformly random in the range (0,0.1). The updating of 
the weights was done by randomly taking one pattern from the training set. For overall 
calculations we used T = 1 and rj = 0.01. 

Since for the training and test data sample both input {x} and correct output {y} 
have to be known for each event, the adjustment of weights and thresholds depends on 
simulated air shower events. For the creation of the training and test showers, we used 
the Monte Carlo code "CORSIKA+ QGSJET" discussed above. The Monte Carlo showers 
were divided into a training sample and test sample and ANN was trained to increase the 
capability for separating the proton-induced events from others. The separation power of 
protons from others may depend upon the chemical composition of primary particles so that 
we trained the ANN using both data samples obtained from the HD and PD primary models 
and checked the difference between them. 

V. RESULTS AND DISCUSSIONS 
A. Behavior of burst events 

First we discuss the behavior of burst events. In Fig. |^, we present the burst size {J2 ^b) 
spectrum observed in our experiment and compare it with the simulation results obtained 
by three different models. 

This figure shows that the CORSIKA+HD and COSMOS+HD models are almost con- 
sistent with the experiment. It is noted that two hadronic interaction models, QGSJET 
in CORSIKA and quasiscaling in COSMOS, can fairly well reproduce many data obtained 
by accelerator and cosmic ray experiments. However, the absolute intensity by the COR- 
SIKA+PD model gives results about three times as high as that by the HD model. This 
difference can be mostly attributed to the difference of the proton flux in both models since 
most selected events are induced by protons, in other words, the observed flux of the burst 
events is very sensitive to the absolute intensity of primary protons. 
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The distribution of the number of fired burst detectors and the air-shower size spectrum 
are also shown in Figs. ^ and ^, respectively, where the experimental results are compared 
with the simulations obtained by the CORSIKA+HD and COSMOS+HD models. From 
these comparisons, we can assure that almost all behavior of the burst events observed are 
compatible with a heavy enriched primary composition at energies around the knee. In 
the previous paper we also discussed the detailed features of the burst events whose 
primary energies are in the knee energy region, say higher than 10^ TeV and reached the 
same conclusion. Based on these results, in the following we try to obtain the primary 
proton spectrum from the observed burst events using the ANN discussed above. 

B. Selection of proton-induced events with ANN 

We trained and tested the ANN using the simulation events obtained from the COR- 
SIKA+HD model, since this model can explain well the behavior of the observed burst 
events as discussed above. For this, 2 x 10'^ events by protons and 2 x 10^ events by other 
nuclei were used as the training data set and the same number of events as the test data set. 
The target value for protons was put to and for other nuclei to 1. A strict middle-point 
condition was used to measure the classification ability of the network, that is, when the 
ANN output is smaller than 0.5, the event is assigned as a proton origin, while when the 
ANN output is larger than 0.5, the event is considered to be an origin of other nuclei. The 
fraction of correct classifications as a function of the number of epochs of the weight updat- 



ing is shown in Fig. |I0[ The dashed and solid lines are for the training and the test data 
sets, respectively. The learning of the network becomes very stable after 300 epochs and the 
change of the weights is small. It is found that the network is able to correctly select 75.7 
% of the two kinds of events we input. The wrong classifications are approximately equally 
distributed among those two. 

As discussed in Sec. Ill, different primary models give different fractions of the events 
produced by each species of primary particles, thus we need to use different values for cutting 
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the network output in order to reduce the wrongly classified events to the desired amount. 
The ANN output distribution of the test events in the case of the HD model is presented 



in Fig. [IT]. It is seen that the proton-induced events can be clearly separated from others with 
a proper cut value of the ANN output. Shown in Fig. |T2| are the ratio of A^(< Vout) /Ntotai 
and the selection efficiency of proton events as a function of the cut Uout in the network 
output, where A^(< yout) is the number of events with the cut < yout and Ntotai is the total 
number of test events used. Here we examined three cases : (1) both training and test data 
sets consist of HD events ; (2) both training and test data sets consist of PD events ; and (3) 
training data set consists of PD events while the test data set consists of HD events. As seen 



in Fig. |T2|, it is confirmed that the ANN training is almost independent upon the primary 
composition and the selection efficiency of proton-induced events is about 90 % when the 
cut value of ANN output yout is set to 0.15. 

Using the ANN trained by the CORSIKA + HD events, we selected 820 candidate events 
induced by protons out of 4274 observed events. 

C. Proton spectrum 

The primary energy of each event can be estimated from a value of Yl, Nb observed with 



the burst detector array. Shown in Fig. 13 is the scatter plots between the burst size 



J2 A"fe and the primary energy Eq of proton-induced events which were selected from the 
data set of the CORSIKA+HD events by setting the ANN output value to 0.15. A fairly 
good correlation between Eq and J2 Nb, as seen in this figure, enables us to estimate the 
primary energies of observed burst events with small ambiguity. A good correlation is also 
found between the air-shower size and the primary energy Eq [ it can be expressed as 
Eq ^ 2.5{GeV) x A^e for A^e > 10^ ], and it is checked that both give almost the same 
values on the primary energy, while the size estimation becomes worse for air-shower events 
with Ne < 10^ because of small number of detectors to be used for fitting. In the present 
analysis, then, we used the burst sizes for the estimation of primary energies. The systematic 
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error on the primary energy estimation was evaluated by the Monte Carlo simulation and 



was estimated to be about 30 % at energies around 500 TeV. Shown in Fig. ^ is the 
effective collecting area of the burst array calculated for primary protons incident at the top 
of the atmosphere isotropically within the zenith angle smaller than 45°. The burst events 
satisfying the selection criteria discussed in Sec. II B and accompanying air showers with 
Np > 10^'^ are selected in this calculation. 



In Fig. |T5|, we present the primary proton spectrum obtained from the burst events, 
which were selected using the ANN trained by the CORSIKA+HD events. To examine 
whether or not the result depends on the primary composition model used, the following 
check was done. For this, first we trained the ANN by using the events obtained from the 
CORSIKA+PD model. Then we selected the proton-induced events from the experimental 
data to obtain the proton spectrum. The primary proton spectrum, thus obtained, is also 



shown in Fig. |l^ to compare with that obtained from the HD composition. Note that in 
spite of a big difference between the HD and PD models on the power index and absolute 
flux of proton component, both results give the same spectrum for protons, as seen in Fig. 
r5| . Hence, we may say that the primary proton spectrum obtained from our experiment 
using the ANN method is almost independent of the primary composition model used in 
the simulation, and it is estimated that the ANN can select the proton-induced events from 
others with an uncertainty of about 10 % under our experimental condition. 

The proton spectrum obtained from this experiment can be represented by the power-law 



fit as shown in Fig. |T5|. The power indexes are estimated to be — 2.97±0.06 and —2.99 ±0.06 
for the spectra obtained using the ANN trained by the CORSIKA+HD and CORSIKA+PD 
events, respectively, where errors quoted are statistical ones. 

It is known that the interpretation of air shower measurements depends on the model 
of the shower development in the atmosphere. The largest uncertainties may originate from 
the hadronic interaction which is not well known at very high energies as well as small 
momentum transfers. Thus, using different hadronic interactions may lead to different pre- 
dictions for some air-shower observables. No drastic changes, however, have been observed 
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on the hadronic interactions at least up to pp collider energies, corresponding to ~ 1000 
TeV in the laboratory system. Also, it is noted that the air shower size observed at high 
altitude weakly depends on the model, while the difference becomes larger near sea level 



Furthermore, we examined in the previous paper that both CORSIKA (QGSJET) and 
COSMOS simulation codes give almost the same results on the behavior of the burst events 
observed with our detector, resulting in that the spectrum obtained here does not depend 
on the simulation code we used. Consequently, we estimate that the systematic errors on 
the proton flux are smaller than 40 % in this experiment. 

Direct measurements of the proton spectrum in the energy region up to about 100 TeV 



T7| , |29|j30| , while statistics is still scanty, may suggest a slightly flat spectrum with the slope 
of -2.5 — -2.7. When both results are combined, we may say that the proton spectrum 
changes its slope at energy around 100 TeV. This may be in favor of shock acceleration 
at SNRs and when we compared this with the all-particle spectrum obtained by the Tibet 



air-shower array |T^, the primary composition becomes heavy dominant at energies around 
the knee. 

D. On the helium spectrum 

Our experiment is also sensitive to the helium component. In order to estimate the 
primary helium spectrum from our experimental data, we adopted the following method. 
Monte Carlo events induced by protons and helium nuclei are flrst gathered as one group and 
its ANN target output is assigned to be 0, while the events induced by other nuclei belong to 
another group with the ANN target output being 1. After training the ANN with the Monte 
Carlo events, then the proton-|-helium events were selected with a proper cut of the ANN 



output as described in our previous paper |21|. The J^^b spectrum of the proton -(-helium 



events minus that of the proton events should give the pure helium spectrum. 

Calculating the effective area for observing the helium-induced events with our burst 
detectors and also using a relation between the burst size J2 and the primary helium 
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energy calculated by the CORSIKA+HD model, we obtained the energy spectrum of primary 



helium nuclei in the energy region above about 100 TeV/n, which is shown in Fig. |T6|. The 
spectrum obtained based on the CORSIKA+PD model is also shown in the same figure to 
compare with each other. Our data is compatible with those extrapolated from the RUNJOB 



5I| and MUBEE |^ data, and the spectrum is not so hard as the JACEE data []T^ at high 



energies. 



VI. SUMMARY 

We have been successfully operating a hybrid experiment of burst detector, emulsion 
chamber and Tibet-II air-shower array since 1996. Using the data obtained with the burst 
detector array and the air-shower array and applying a neural network analysis to this data 
set, we obtained the energy spectrum of primary protons in the energy range from 200 to 
1000 TeV. The spectral index is estimated to be —2.97 ± 0.06, suggesting that the proton 
spectrum should steepen at energies of 100 TeV when compared with direct observations 
done in the lower energy region. 

We also estimated the primary helium spectrum at particle energies around 1000 TeV, 
which may have almost same spectral slope with the proton spectrum, though the statistics 
is still not enough. 

Using gamma family events, those observed with the emulsion chamber, accompanied 
by air showers, we can estimate the primary proton spectrum in the energy region from 
10^ TeV to ~ 10'^ TeV and the result will be reported in very near future [^. Then, the 
Tibet air-shower experiment can measure the primary proton spectrum in the wide energy 
range from 200 TeV to ~ 10^ TeV and provide vital information necessary for clarifying the 
acceleration mechanism of cosmic ray particles at very high energies. 
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APPENDIX A: PERFORMANCE OF THE BURST DETECTOR 

Each burst detector contains a plastic scintillator with the size of 160cm x 50cm x 2cm. 
A PIN PD (HPK S2744-03) with an effective area of 2cm x 1cm was equipped at each of 



four corners of the scintillator, as shown in Fig. To detect signals from a PD for burst 
particles ranging from 10^ to 10^, a preamplifier with an amplification factor of 260 operating 
in the frequency range from ITkHz to 44 MHz (current- current type) was developed. An 
ADC value from each PD, depending on the size and the hit position of a burst (shower) 
fallen in the burst detector, can be expressed as KNb{r), where r is the distance between 
a PD and the burst position in the scintillator, Ni, is the burst size, and i^' is a constant. 
Using the ADC values from four corners, we can estimate the size and hit position for each 
burst event using a least-squares method. In this formula, /(r) denotes the attenuation of 
photons in the scintillator. In general /(r) can be expressed as exp(-r/A) except at small 
distance r and A takes a value around 350 cm for the present scintillator. Since the size of 
the burst detector is smaller than the attenuation length, errors of the estimation of burst 
hit position become very large. So we first slightly polished one face of each scintillator with 
rough sandpaper (No. 60) to make photons scatter randomly on this face. Then we found 
that /(r) can be well approximated as and a ~ 1.1 ~ 1.2. This relation was confirmed 
by using a nitrogen gas laser and also cosmic ray muons. This dependence on the distance 
r is sufficient to estimate the burst position in the detector. 
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We also installed a calibration unit which consists of four blue light-emitting diodes 
(LED's) each having a peak wave length of 450 nm. The LED unit is put on the center of 
each scintillator and is illuminated to transmit light through the scintillator to each PD at 
the corner uniformly, and then all the ADC's are calibrated at every 10 min for actual run. 
This calibration system provides information about a relative change of ADC values, which 
may cause a large error for the estimation of burst hit positions and burst sizes. 

We examined the performance of the burst detector using electron beams of 1.0 GeV/c 
from the KEK-Tanashi Electron Synchrotron. The electron beams, ranging from several 
xlO^ to ~ 3 X 10^ per pulse, were vertically exposed to various positions on the surface of 
the burst detector. 



Figure |T8| shows the dependence of the ADC values on the distance r, obtained with the 
electron beams, where r is the distance between the beam hit position and PD. The result 
can be well fitted by a power law of r, where the number of incident electrons measured by 
the probe scintillator was normalized to 10^ particles. 

Using the ADC values from four PD's, the beam positions exposed on the face of the 
detector and its intensities (number of electrons) were estimated to compare with the true 
ones. The distribution of the difference between estimated and actual beam positions is 



shown in Fig. |I9|. We present scatter plots of the estimated number and irradiated number 



of electrons in Fig. |20, and the distribution of the ratio between them is shown in Fig. 



2T| . From these figures, it is concluded that the hit position of a burst in each detector can 
be estimated with an inaccuracy of less than 10 cm and errors for the size estimation are 
smaller than 10 % for the bursts with size > 10^ particles. 

APPENDIX B: PRIMARY COSMIC RAY COMPOSITION 



The energy spectra of respective components assumed in the heavy dominant (HD) 
model and proton dominant (PD) models are shown in Figures ^ (a) and (b), respectively. 
The all-particle spectra obtained by the experiments : Tibet [|T7|, PROTON satellite 
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JACEE ||3J] , and AKENO [^] are plotted in both figures. The all-particle spectrum in each 
model is normalized to the Tibet data at energies around the knee. 
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FIG. 1. Arrangement of 100 burst detectors set up in two rooms. The area of each burst 
detector is 50 cm x 160 cm and four emulsion chambers are set up on each burst detector. 
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FIG. 2. Schematic side view of each unit of emulsion chamber and burst detector. High 
sensitive x-ray films are inserted at every 2 c.u. in emulsion chamber. Total thickness of lead 
plates is 15 c.u. (7.5 cm). 
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FIG. 3. Example of air-shower core event observed in the burst detectors. Rhombi denote the 
size of events observed in each burst detector and its area is logarithmically proportional to the 
burst size. 




2.0 2.5 3.0 3.5 4.0 4.5 5.0 
log1 0(burst size in OTHER-section) 



FIG. 4. Burst size distribution in the OTHER section. The five curves denote the different 
positions of the TOP detector being in the first, second, third, fourth, and fifth column in the 
TOP section, respectively. 
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FIG. 5. Scatter dots of the total burst size J2 and the shower size Ne for the CORSIKA+HD 
simulation events induced by protons (a) and other nuclei (b). 



28 




FIG. 6. Feed-forward neural network with one hidden layer. 



10-^ 



r-i 
< 

10-3 

I 10-^" 
< 



5.0 5.5 6.0 6.5 7.0 

LoglOlZBurst Size) 

FIG. 7. Burst size spectrum. The open circles, solid line, long-dashed line, and dotted line 
denote the experimental data, CORSIKA+HD, CORSIKA+PD, and COSMOS+HD simulation 
events, respectively. 
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FIG. 8. Number distribution of the fired burst detectors for each event. The open circles, 
histogram, and dotted line denote the experimental data, CORSIKA+HD and COSMOS+HD 
simulation events, respectively. 
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FIG. 9. Size distribution of the air showers accompanied by the burst events with Y,^b> 10^- 
The open circles and solid line denote the experimental data and CORSIKA+HD simulation events, 
respectively. 
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FIG. 10. Network performance as a function of the number of training epochs, where the dashed 
and soUd lines are for the training data set and test data set by HD model, respectively. 
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FIG. 11. ANN output distribution of the CORSIKA+HD simulation events. The solid and 
dotted lines denote the events induced by protons and other nuclei, respectively. 
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FIG. 12. Ratio of the number of selected events with the network out smaller than yout to the 
total test events and the selection efficiency of the proton-induced events, expressed as a function 
of the network output yout- Solid line, dashed line, and dotted line are the cases where the training 
set is the HD data and the test set is the HD data, where the training set is the PD data and the 
test set is the PD data, and where the training set is the PD data and the test set is the HD data, 
respectively. 
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FIG. 13. Scatter plots between the primary energy Eq versus the total burst size X] for 
proton-induced events in the case of the CORSIKA+HD model. The events are selected by setting 
the ANN output value to 0.15. The solid circles denote the average values, and the solid line is a 
fit by the relation £^0 = 1200{J2 Nb/lO^ f-^^ (TeV). 
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FIG. 14. Effective collecting area of the burst array for primary protons entering isotropically 
at the top of atmosphere (zenith angle < 45°). For the selection criteria of the burst events, see 
text. 



33 



10^ 



(N 

< 

> 

a 10^ 



* 

(N 
< 

o 



10^ 



O JACEE 
□ RUNJOB 
A MUBEE 



• This work by CORSIKA+HD 
■ This work by CORSIECA+PD 



Our data fit, index=-2.97 



10 100 1000 

Primary Proton Energy (TeV) 



FIG. 15. Energy spectrum of primary protons. The filled circles and squares stand for the 
experimental results obtained using the ANN's trained by the CORSIKA+HD and CORSIKA+PD 
events, respectively. Our results are compared with other direct measurements by JACEE [10], 
RUNJOB [30], and MUBEE [31]. The dashed line is a best fit to our data. 
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FIG. 16. Energy spectrum of primary helium nuclei. For details, see text. Our results are 
compared with other direct measurements by JACEE [10], RUNJOB [30], and MUBEE [31]. 
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FIG. 17. Schematic view of the burst detector used in this experiment. Numerals shown in 
the figure are in units of mm. 



10 ^ I I I 

■ ■■ exp.dat 1^9—1 

0.563*(x/83.8)**-1.12 



10^ r ' 




■^Q-^l I I I I 

10° 10^ 10^ 10^ lo'' 

Distance (cm) 



FIG. 18. Attenuation of photons in the scintillator used for the burst detector, obtained using 
electron beams. 
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FIG. 19. Distribution (integral) of the difference between estimated and irradiated positions. 
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FIG. 20. Scatter plots of estimated and irradiated number of electrons. The number of 
electrons at various beam positions on the face of the detector is normalized to 10^ electrons. 
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FIG. 21. Distribution of the ratio of estimated and irradiated number of electrons shown in 



Figure 21. Dotted hne is a Gaussian fit. 
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FIG. 22. Primary cosmic ray composition for (a) the HD model and (b) the PD model. The 
all-particle spectrum, which is a sum of each component, is normalized to the Tibet data. 
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